Random close packing of polydisperse hard spheres 
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We study jammed configurations of hard spheres as a function of compression speed using an 
event-driven molecular dynamics algorithm. We find that during the compression, the pressure 
follows closely the metastable liquid branch until the system gets arrested into a glass state as the 
relaxation time exceeds the compression speed. Further compression yields a jammed configuration 
that can be regarded as the infinite pressure configuration of that glass state. Consequently, we find 
that the density of jammed packings varies from 0.638 to 0.658 for polydisperse hard spheres and 
from 0.635 to 0.645 for pure hard spheres upon decreasing the compression rate. This demonstrates 
that the density at which the systems falls out of equilibrium determines the density at which the 
system jams at infinite pressure. In addition, we give accurate data for the jamming density as a 
function of compression rate and size polydispersity. 

PACS numbers: 82.70.Dd,61.43.Fs,64.70.Q- 



Dense packings of hard particles are relevant for many 
applications ranging from the processing of granular ma- 
terials to the development of new materials. In addi- 
tion, hard spheres have provided a good starting point in 
the study of liquids, glasses, crystals, colloids, and pow- 
ders. The study of dense packings of hard spheres has a 
long history, dating back to Kepler for crystalline pack- 
ings and Bernal 1] for random packings. In 2005, Hales 
proved Kepler's conjecture 0] that the densest packing 
of identical hard spheres is achieved by the stacking of 
close-packed hexagonal planes yielding a packing frac- 
tion </> = ira 3 N/6V rj 0.74. However, when a system of 
hard spheres is compressed quickly it does not reach this 
maximum density and it jams at lower densities. Jam- 
ming phenomena are generic since atomic, colloidal and 
granular systems can all be jammed in a state far out of 
equilibrium by quickly cooling, compressing or unloading 

Many authors speculate that the packing fraction of 
the densest amorphous hard-sphere configuration is well- 
defined 0, [B| although its precise value is not known as 
random-close-packing densities 4>rcp ranging from 0.634 
[H to 0.648 0] are reported. On the other hand, Torquato 
et al argued that ^ r cp is ill-defined as the density 
can always be increased by introducing crystalline order 
into the system. Very recently, it is speculated and/or 
deduced from mean-field models that jammed configu- 
rations can be regarded as the infinite pressure limit of 
glassy states [H, [t| . The reason is that the glass phase for 
each density on the metastable fluid branch consists of 
a certain group of glassy states, which will all follow the 
same metastable glass branch up to the same jamming 
density by fast compression @, Q . 

In this letter, we study jammed configurations of hard 
spheres as a function of compression speed. We show 
that denser amorphous packings can be generated even 
without introducing any crystalline order by lowering 
the compression speed. However, in the case of iden- 



tical spheres, the system (partially) crystallizes for suffi- 
ciently slow compressions. Crystallization can be avoided 
by introducing size polydispersity. A small size poly- 
dispersity affects 0rcp only slightly but hampers crystal 
growth severely. This allows us to study jamming for a 
wider range of compression rates. We find volume frac- 
tions ranging from 0.638 till 0.658 for hard spheres with 
10% size polydispersity upon decreasing the compression 
rates. Additionally, we give accurate data for the jam- 
ming density (j)j as a function of compression rate and 
size polydispersity, which is important for experiments 
on colloidal systems. The packing fraction of a poly- 
disperse colloidal system is often determined by setting 
the packing fraction of a centrifuged sedimen t II 011 equal 
to the i^rcp obtained from simulations ll|, ll2j . How- 
ever, these simulation results are very inaccurate and do 
not take into account any compression rate dependence, 
casting doubts on the accuracy of the volume fractions 
determined in experiments via this route. As our re- 
sults show that the jamming density strongly depends on 
the compression rate, we prefer to denote the density of 
the jammed configuration with the jamming density </>j 
rather than with the random-close-packing density <^rcp- 

We perform simulations with the event driven molec- 
ular dynamics algorithm of Ref. [l3(. Modifications as 
described by Speedy Q were used to fix the tempera- 
ture of the system and to define a compression speed 
r = da/dt, where we use the MD time as our unit of 
time. In addition the algorithm was adapted to keep the 
polydispersity constant during the particle growth 14 1 . 
The polydispersity was sampled from a log-normal dis- 
tribution. The log-normal distribution is nearly identical 
to a normal (Gaussian) distribution for small polydisper- 
sities but has the advantage that it is zero for negative 
diameters. 

We determine the equation of state (EOS) from simu- 
lations of 2000 particles with a size polydispersity of 10%. 
We average our results over 50 different runs. To check 
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FIG. 1: a) Pressure [3Pa 3 as a function of packing fraction 
<f> for a system of hard spheres with 10% size polydispersity 
and for varying compression rates F as labeled. The red and 
blue (dark and light full) line denote the equilibrium equation 
of state for the fluid and solid phase, respectively, while the 
horizontal line denote fluid-solid coexistence. The dotted lines 
denote the fits to the simulation data using Eq. (1). b) The 
same results as in a) but now plotted as a function of cf> and 
the inverse pressure l/(f3Pa ) so that the infinite pressure 
limit is clearly visible. 



for finite size effects we perform simulations with up to 
200000 particles and we find good agreement within the 
statistical accuracy. We plot the EOS for varying com- 
pression rates in Fig. [1^ along with the equilibrium EOS 
for the fluid [IH and the solid phase [l6| . We checked for 
crystallization during our simulations, but we did not find 
any crystalline order in any of our compression runs. We 
observe that the pressure follows initially the equilibrium 
EOS of the fluid phase until the system becomes arrested 
into a glass state as the relaxation time of the system ex- 
ceeds the compression rate. At this density the pressure 
increases much faster than that of the equilibrium fluid 
EOS upon further compression. We find that the density 
of the jammed configuration at infinite pressure increases 
with slower compressions. The reason is that the system 
has more time to equilibrate for slow compressions, and 
hence the system falls out of equilibrium at a higher den- 



sity on the metastable fluid branch. Further compres- 
sion of this glass phase yields a jammed configuration 
with a density that is higher than for fast compressions. 
We find clear evidence that the jamming density of the 
amorphous packing depends strongly on the glass transi- 
tion density, i.e., the density where the system becomes 
arrested. Therefore we find a finite range of densities 
for the jamming density depending on the compression 
procedure and history of the sample. 

This behaviour is well known for molecular glasses [l7j . 
Indeed, Fig. QJi resembles the picture that is found for 
molecular glasses with <f> and the inverse of (3Pa 3 play- 
ing the role of the inverse of the specific volume and the 
temperature, respectively. To this end, we plot our re- 
sults in Fig. [TJd in the <f> - l/((3Pa 3 ) representation, 
where (3 — l/kgT with ks Boltzmann's constant. Wc 
now find striking similarities with the sketched phase di- 
agrams in Refs. [1, @|. Recent theoretical calculation 
using the replica method by Parisi and Zamponi predict 
indeed that different glasses can jam at different densi- 
ties upon compression, and that the pressure of the glass 
phase close to jamming is well-described by a power law 
(3P/ p oc 1/{4>j — 4>) with 4>j the jamming density at infi- 
nite pressure [8|. We observe in Fig. [TJd an almost linear 
behavior for the inverse pressure as a function of <f> for 
the glass phase, which we can fit remarkably well over 
the full range using the free volume scaling [j]: 



6 1/3 6 2 ' 3 
(3Pa 3 = a ^ J 



(i) 



where a and <f>j are fitting parameters. The leading or- 
der term of Eq. (1) yields the power law as predicted in 
[8], close to jamming. In addition, the theory predicts 
an ideal glass transition at a Kauzmann packing frac- 
tion (j)K — 0.617 on the metastable fluid branch, yielding 
a jammed configuration upon compression of this ideal 
glass with a random close packing density of 4> rcp = 0.683 
[8j. However, the existence of a thermodynamic ideal 
glass transition is still heavily debated. We note that our 
results do not depend on the existence or non-existence 
of such a glass transition. Our results show a nonequilib- 
rium glass transition at densities in the volume fraction 
range of 0.50-0.59, which is far below the theoretical pre- 
dictions for the ideal glass. This is to be expected as 
the structural relaxation time diverges on approaching 
the ideal glass transition Hence, it is impossible to 
reach <j>K, since already at lower densities, the fluid gets 
arrested in a noncquilbrium glass as the relaxation time 
becomes longer than the simulation time. 

In order to investigate the divergence of the structural 
relaxation time on approaching (f>x, and to estimate what 
the maximum density is at which we can still equilibrate 
a state on the metastable fluid branch, we perform con- 
stant volume simulations of a system of 20000 spheres 
with a 10% size polydispersity and varying packing frac- 



3 



100 



o. 



0.01 



0.001 




1 10 100 1000 10000 100000 

t 

FIG. 2: The pressure difference /3(P — Pcs)cr 3 with respect to 
the Carnahan-Starling equation of state Pes obtained from 
constant volume simulations as a function of MD time for a 
system of hard spheres with 10 % size polydispersity. The 
runs are started with an initial configuration obtained from a 
fast compression. 



tions. We start the runs with an initial configuration ob- 
tained from a fast compression. Fig. [2]shows the pressure 
difference with respect to the Carnahan-Starling EOS as 
a function of time. We clearly see that the pressure ini- 
tially decays towards an intermediate state for all <f>. Sub- 
sequently, large collective rearrangements are required to 
relax the system further. Finally, we observe that the 
pressure reaches the value predicted by the equilibrium 
Carnahan-Starling EOS. The time scale for the system to 
equilibrate to the equilibrium fluid phase is comparable 
to the relaxation time r a that can be determined from 
the decay of the intermediate scattering function [l^ |. 
We clearly observe in Fig. [5] that the relaxation time di- 
verges with packing fraction. The equilibration time of 
a system with 20000 particles at <f> = 0.585, costed more 
than 10 5 MD time steps, which is equivalent to two weeks 
on a desktop PC. The equilibration time for cf> — 0.59 is 
expected to be more than 20 weeks. Hence, the ideal 
glass transition at 4>k, and the corresponding random 
close packing that can be achieved by compressing the 
ideal glass, are both inaccessible. Instead the system will 
fall out of equilibrium into a non-equilibrium glass state 
at a density that depends strongly on the compression. 

We therefor study systematically the jamming density 
as a function of compression speed for several size poly- 
dispersities 6. We perform simulations of 2000 particles 
using varying compression rates and we terminate the 
simulations when the time between successive collisions 
becomes of the same order of magnitude as our numer- 
ical accuracy, yielding pressures (3 Pa 3 up to 100000 for 
slow compressions. To determine the jamming density 
<j)j, we fit the equation of state close to jamming (the 
last few volume percent) with Eq. (1). We average our 
results over 50 different runs. Fig. [3^ shows the jamming 
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FIG. 3: a) Jamming density <f>j as a function of the inverse 
compression rate T" 1 for polydispersities 5 as reported in the 
labels. The lines denote fits to our results (see text). For 
monodisperse hard spheres, the line ends where (partial) crys- 
tallization is observed, b) Jamming density </>j for = 20 
and 10 4 as a function of size polydispersity S. The data from 
Nolan et al. [ll|] and Schaertl et al. [l^] are denoted by the 
crosses and plusses, respectively. 



density as a function of the inverse compression rate for 
varying polydispersities S. We have fitted these extrap- 
olations with (f> = a + 6/log(r), where a and b are fit- 
ting parameters. The jamming density <pj for pure hard 
spheres ranges from 0.635 - 0.645 for 10 < T" 1 < 1000. 
For faster compressions < 10, the simulations do not 
yield jammed configurations and the system (partially) 
crystallizes for > 1000 as can be observed in Fig. f3]as 
<f>j increases rapidly. We also observe that the jamming 
density 4>J increases with increasing size polydispersity. 
For S = 10%, we find that cj>j varies from 0.638 - 0.658. 
As the size polydispersity prevents crystallization, we are 
now able to study <j)j for compression rates across five 
orders of magnitude. Although the slope of the curves 
decreases with increasing T _1 (slower compressions), it 
is hard to justify an extrapolation to infinitely slow com- 
pression rates. Figure Wp shows the jamming density as 
a function of size polydispersity 8 for a compression rate 
r _1 = 20 and 10 4 . For comparison, we also plot the data 
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from Nolan 11 1 and Schaertl Our results are close 



to those of Nolan [111 ], The results from Schaertl et al. 



deviate from our results but can be explained since they 
have done single runs with low accuracy fl2| . The strong 
dependence of the jamming density on the compression 
speed explains the range of densities that has been found 
in the literature, obtained by different authors and by 
using different algorithms. Our results also show that a 
size polydispersity of up to 5% does not increase the jam- 
ming density significantly from the monodisperse case. 
A much larger dependence on the polydispersity is often 
used in the experiments [13] based on the results of Ref. 
[13], casting doubts on the precise values for the volume 
fractions determined in experiments via this route. 

In conclusion, we have studied the jamming density 
of hard spheres as a function of compression speed. We 
find that during the compression, the pressure follows 
closely the metastable liquid branch until the system gets 
arrested into a glass state. Further compression of the 
glass phase yields an almost linear behavior for the in- 
verse pressure as predicted theoretically 0, @], provid- 
ing evidence that the jammed configuration can be re- 
garded as the infinite pressure configuration of that glass 
state. We find, as expected, higher jamming densities for 
slower compression. The reason for this is that the sys- 
tem has more time to equilibrate for slower compressions, 
and hence the density increases at which the system falls 
out of equilibrium. Consequently, we find that the den- 
sity of jammed packings varies from 0.635 to 0.645 for 
pure hard spheres upon decreasing the compression rate. 
For slower compressions we observe partial crystalliza- 
tion and for faster compression our simulations do not 
result in jammed states. For slightly polydisperse (10%) 
hard spheres, the jamming density varies from 0.638 to 
0.658. Additionally, we give accurate data for the jam- 
ming density as a function of compression rate and size 
polydispersity. Our results provide evidence that there 



is a fundamental difference between the ideal glass tran- 
sition Q , which takes place at a fixed density (if it does 
exists at all), and the jamming transition [13], that takes 
place at a range of densities, as was already speculated 
by [1, 13] • However, there is a link between the density at 
which the systems falls out of equilibrium and the density 
at which the system jams at infinite pressure. 
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